#ifndef _MESH_
#define _MESH_

#include "Domain.h"
#include "Dofs.h"
#include "Grid.h"
#include <map>
#include <functional>
#include <Eigen/Sparse>
#include <fstream>
using namespace std;
/**
 * @brief 定义域划分成网格Mesh,Mesh包含小网格Grid,每个Grid的类型由TemplateElement确定
 * 
 * @tparam int DIM 维度
 */
template <unsigned int DIM>
class Mesh
{
public:
    /**
     * @brief 默认构造函数
     */
    Mesh() = default;

    /**
     * @brief 默认析构函数
     */
    virtual ~Mesh() = default;
    //virtual Mesh(Domain<DIM>,vector<int>) = 0;
    /**
     * @brief 返回自由度的总个数
     * @details _GridList.size()*getNumDofs();
     * @return 
     */
    virtual int getTotalNumDofs() = 0;
    /**
     * @brief 返回网格的总个数
     * @details _GridList.size()
     * @return [description]
     */
    virtual int getNumGrid() = 0;
    /**
     * @brief 返回存储Grid的vector
     * @details _GridList
     * @return [description]
     */
    virtual vector<Grid<DIM>>& getGridVector() = 0;
    /**
     * @brief 返回序号为_index的Grid
     * 
     * @param _index Grid的序号
     * @return _GridList[i];
     */
    virtual Grid<DIM> &getGrid(int _index) = 0;
    /**
     * @brief   返回vector，存储边界上的自由度
     * @details 对所有grid,依据seg找出边界的grid,再确定对应边界的自由度下标是哪几个
     *
     * @param component 试图对边界分段，并标号
     * @return 存储边界上自由度的vector
     */
    virtual vector<Dofs<DIM>>& getBoundaryDofs(int component = 0) = 0;

    virtual void distributeDofs() = 0;
    
    virtual vector<vector<Dofs<DIM>>>& getDofsOnGrid() = 0;

    virtual vector<int>& getSegment() = 0;
protected:
    vector<Grid<DIM>> gridList;/*< 存放所有的Grid*/
    vector<vector<Dofs<DIM>>> dofsOnGrid;/*< 存放grid上的所有自由度 */
    vector<Dofs<DIM>> boundaryDofs;/*< 存放所有边界上的自由度*/
    vector<int> seg;/*< 存放各维度分段数*/
};

class Q1Mesh : public Mesh<2>
{
public:
    Q1Mesh() = default;
    /**
     * @brief 含参构造函数,根据给定的定义域和划分段数来划分网格。
     * @param domain 定义域
     * @param seg 每个维度要划分的段数
     */
    Q1Mesh(Domain<2> domain, vector<int> _seg);
    ~Q1Mesh() = default;
   
    int getTotalNumDofs();
    int getNumGrid();
    vector<Grid<2>>& getGridVector();
    Grid<2>& getGrid(int _index);
    vector<Dofs<2>>& getBoundaryDofs(int component = 0);
    /**
     * @brief 得到分段
     * 
     */
    vector<int>& getSegment();
    /**
     * @brief 给定模板单元的阶数，在每个Grid上划定自由度。
     * @details 仅实现二维一阶、二阶模板单元
     * 
     * @param _order_of_element 模板单元的阶数
     */
    void distributeDofs();
    vector<vector<Dofs<2>>>& getDofsOnGrid();

protected:
    //vector<int> seg;/*< 存放各维度分段数*/
};


Q1Mesh::Q1Mesh(Domain<2> domain, vector<int> _seg)
{
    assert(_seg.size() == 2);
    Point<2> p0 = domain.getVertex(0);
    Point<2> p2 = domain.getVertex(2);
    double xLen = (p2[0] - p0[0])/(1.0*_seg[0]);
    double yLen = (p2[1] - p0[1])/(1.0*_seg[1]);
    for(int j = 0; j < _seg[1]; ++j)
    {
	for(int i = 0; i < _seg[0]; ++i)
	{
	    vector<Point<2>> pnts;
	    pnts.push_back({p0[0]+i*xLen,p0[1]+j*yLen});
	    pnts.push_back({p0[0]+(i+1)*xLen,p0[1]+j*yLen});
	    pnts.push_back({p0[0]+(i+1)*xLen,p0[1]+(j+1)*yLen});
	    pnts.push_back({p0[0]+i*xLen,p0[1]+(j+1)*yLen});
	    Grid<2> grid(pnts);
	    Mesh<2>::gridList.push_back(grid);
	}
    }
    seg = _seg;
    
};


int Q1Mesh::getNumGrid()
{
    return Mesh<2>::gridList.size();
};


Grid<2>& Q1Mesh::getGrid(int _index)
{
    return Mesh<2>::gridList[_index];
};


vector<Grid<2>>& Q1Mesh::getGridVector(){
    return Mesh<2>::gridList;
}


vector<int>& Q1Mesh::getSegment()
{
    return seg;
};


vector<vector<Dofs<2>>>& Q1Mesh::getDofsOnGrid(){
    return Mesh<2>::dofsOnGrid;
}


vector<Dofs<2>>&  Q1Mesh::getBoundaryDofs(int component)
{
    //对所有grid,依据seg找出边界的grid,再确定对应边界的自由度下标是哪几个

    //只实现维数为2 精度为1、2
    //assert(DIM ==2);
    
    for(int k = 0; k < seg[0]; k++)
	Mesh<2>::boundaryDofs.push_back(Mesh<2>::dofsOnGrid[k][0]);
    for(int k = 0; k <= seg[0] * (seg[1] - 1); k = k + seg[0])
	Mesh<2>::boundaryDofs.push_back(Mesh<2>::dofsOnGrid[k][3]);
    for(int k = seg[0] - 1; k < seg[0] * seg[1]; k = k + seg[0])
	Mesh<2>::boundaryDofs.push_back(Mesh<2>::dofsOnGrid[k][1]);
    for(int k = seg[0] * (seg[1] - 1); k < seg[0] * seg[1]; k++)
	Mesh<2>::boundaryDofs.push_back(Mesh<2>::dofsOnGrid[k][2]);
	
    // else if(_order_of_element == 2)
    // {
    // 	for(int k = 0; k < seg[0]; k++){
    //         Mesh<DIM>::boundaryDofs.push_back(Mesh<DIM>::dofsOnGrid[k][0]);
    //         Mesh<DIM>::boundaryDofs.push_back(Mesh<DIM>::dofsOnGrid[k][1]);}
    // 	for(int k = 0; k <= seg[0] * (seg[1] - 1); k = k + seg[0]){
    //         Mesh<DIM>::boundaryDofs.push_back(Mesh<DIM>::dofsOnGrid[k][7]);
    //         Mesh<DIM>::boundaryDofs.push_back(Mesh<DIM>::dofsOnGrid[k][6]);}
    // 	for(int k = seg[0] - 1; k < seg[0] * seg[1]; k = k + seg[0]){
    //         Mesh<DIM>::boundaryDofs.push_back(Mesh<DIM>::dofsOnGrid[k][2]);
    //         Mesh<DIM>::boundaryDofs.push_back(Mesh<DIM>::dofsOnGrid[k][3]);}
    // 	for(int k = seg[0] * (seg[1] - 1); k < seg[0] * seg[1]; k++){
    //         Mesh<DIM>::boundaryDofs.push_back(Mesh<DIM>::dofsOnGrid[k][5]);
    //         Mesh<DIM>::boundaryDofs.push_back(Mesh<DIM>::dofsOnGrid[k][4]);}
    // }
    
    return Mesh<2>::boundaryDofs;
};

//假定所有单元类型、阶数都相同

void Q1Mesh::distributeDofs()
{
    //以下示范DIM=2，_order_of_element = 1 | 2;
    //assert(DIM ==2);
    int k = 0;
    vector<Dofs<2>> _dofs;

    for (Grid<2> &grid : Mesh<2>::gridList)
    {
	    
	_dofs.clear();
	vector<Point<2>> p = grid.GetVertices();
	_dofs.push_back({{p[0][0], p[0][1]}, k + k / seg[0]});
	_dofs.push_back({{p[1][0], p[1][1]}, k + k / seg[0] + 1});
	_dofs.push_back({{p[2][0], p[2][1]}, k + k / seg[0] + seg[0] + 2});
	_dofs.push_back({{p[3][0], p[3][1]}, k + k / seg[0] + seg[0] + 1});
	Mesh<2>::dofsOnGrid.push_back(_dofs);
	k++;	    
    }

    // else if(_order_of_element == 2)
    // {
    // 	for (Grid<DIM> &grid : Mesh<DIM>::gridList)
    //     {
    //         _dofs.clear();
    //         int a = seg[0];
    //         vector<Point<DIM>> p = grid.GetVertices();
    //         _dofs.push_back({{p[0][0], p[0][1]}, k / a * (4 * a + 2) + k % a * 2});
    //         _dofs.push_back({{(p[0][0] + p[1][0])/2.0, (p[0][1] + p[1][1])/2.0}, k / a * (4 * a + 2) + k % a * 2 + 1});
    //         _dofs.push_back({{p[1][0], p[1][1]}, k / a * (4 * a + 2) + k % a * 2 + 2});
    //         _dofs.push_back({{(p[1][0]+p[2][0])/2.0,(p[1][1] + p[2][1])/2.0}, k / a * (4 * a + 2) + 2 * a + 1 + k % a * 2 + 2});
    //         _dofs.push_back({{p[2][0],p[2][1]}, (k / a + 1) * (4 * a + 2) + k % a * 2 + 2});
    //         _dofs.push_back({{(p[2][0]+p[3][0])/2.0,(p[2][1]+p[3][1])/2.0}, (k / a + 1) * (4 * a + 2) + k % a * 2 + 1});
    //         _dofs.push_back({{p[3][0],p[3][1]}, (k / a + 1) * (4 * a + 2) + k % a * 2});
    //         _dofs.push_back({{(p[3][0]+p[0][0])/2.0,(p[3][1]+p[0][1])/2.0}, k / a * (4 * a + 2) + 2 * a + 1 + k % a * 2});
    //         _dofs.push_back({{(p[0][0]+p[2][0])/2.0,(p[0][1]+p[2][1])/2.0}, k / a * (4 * a + 2) + 2 * a + 1 + k % a * 2 + 1});
    // 	    Mesh<DIM>::dofsOnGrid.push_back(_dofs);
    //         k++;
    //     }
    // }
    // else
    // {
    //     std::cout <<"抱歉，干不了。\n"; 
    // }
};

int Q1Mesh::getTotalNumDofs()
{
    //此为Q_1_E算法，Q_2_E算法
    return (seg[0] + 1) * (seg[1] + 1);
    // else if(_order_of_element ==2)
    // 	return (2 * seg[0] + 1) * (2 * seg[1] + 1);
    // else
    // 	return 0;
};

class Q2Mesh : public Mesh<2>
{
public:
    Q2Mesh() = default;
    /**
     * @brief 含参构造函数,根据给定的定义域和划分段数来划分网格。
     * @param domain 定义域
     * @param seg 每个维度要划分的段数
     */
    Q2Mesh(Domain<2> domain, vector<int> _seg);
    ~Q2Mesh() = default;
   
    int getTotalNumDofs();
    int getNumGrid();
    vector<Grid<2>>& getGridVector();
    Grid<2>& getGrid(int _index);
    vector<Dofs<2>>& getBoundaryDofs(int component = 0);
    /**
     * @brief 得到分段
     * 
     */
    vector<int>& getSegment();
    /**
     * @brief 给定模板单元的阶数，在每个Grid上划定自由度。
     * @details 仅实现二维一阶、二阶模板单元
     * 
     * @param _order_of_element 模板单元的阶数
     */
    void distributeDofs();
    vector<vector<Dofs<2>>>& getDofsOnGrid();

protected:
    //vector<int> seg;/*< 存放各维度分段数*/
};


Q2Mesh::Q2Mesh(Domain<2> domain, vector<int> _seg)
{
    assert(_seg.size() == 2);
    Point<2> p0 = domain.getVertex(0);
    Point<2> p2 = domain.getVertex(2);
    double xLen = (p2[0] - p0[0])/(1.0*_seg[0]);
    double yLen = (p2[1] - p0[1])/(1.0*_seg[1]);
    for(int j = 0; j < _seg[1]; ++j)
    {
	for(int i = 0; i < _seg[0]; ++i)
	{
	    vector<Point<2>> pnts;
	    pnts.push_back({p0[0]+i*xLen,p0[1]+j*yLen});
	    pnts.push_back({p0[0]+(i+1)*xLen,p0[1]+j*yLen});
	    pnts.push_back({p0[0]+(i+1)*xLen,p0[1]+(j+1)*yLen});
	    pnts.push_back({p0[0]+i*xLen,p0[1]+(j+1)*yLen});
	    Grid<2> grid(pnts);
	    Mesh<2>::gridList.push_back(grid);
	}
    }
    seg = _seg;
    
};


int Q2Mesh::getNumGrid()
{
    return Mesh<2>::gridList.size();
};


Grid<2>& Q2Mesh::getGrid(int _index)
{
    return Mesh<2>::gridList[_index];
};


vector<Grid<2>>& Q2Mesh::getGridVector(){
    return Mesh<2>::gridList;
}


vector<int>& Q2Mesh::getSegment()
{
    return seg;
};


vector<vector<Dofs<2>>>& Q2Mesh::getDofsOnGrid(){
    return Mesh<2>::dofsOnGrid;
}


vector<Dofs<2>>&  Q2Mesh::getBoundaryDofs(int component)
{
    //对所有grid,依据seg找出边界的grid,再确定对应边界的自由度下标是哪几个

    //只实现维数为2 精度为1、2
    //assert(DIM ==2);
    
    // for(int k = 0; k < seg[0]; k++)
	// Mesh<2>::boundaryDofs.push_back(Mesh<2>::dofsOnGrid[k][0]);
    // for(int k = 0; k <= seg[0] * (seg[1] - 1); k = k + seg[0])
	// Mesh<2>::boundaryDofs.push_back(Mesh<2>::dofsOnGrid[k][3]);
    // for(int k = seg[0] - 1; k < seg[0] * seg[1]; k = k + seg[0])
	// Mesh<2>::boundaryDofs.push_back(Mesh<2>::dofsOnGrid[k][1]);
    // for(int k = seg[0] * (seg[1] - 1); k < seg[0] * seg[1]; k++)
	// Mesh<2>::boundaryDofs.push_back(Mesh<2>::dofsOnGrid[k][2]);
            
            
                    for(int k = 0; k < seg[0]; k++){
                        Mesh<2>::boundaryDofs.push_back(Mesh<2>::dofsOnGrid[k][0]);
                        Mesh<2>::boundaryDofs.push_back(Mesh<2>::dofsOnGrid[k][1]);}
                    for(int k = 0; k <= seg[0] * (seg[1] - 1); k = k + seg[0]){
                        Mesh<2>::boundaryDofs.push_back(Mesh<2>::dofsOnGrid[k][7]);
                        Mesh<2>::boundaryDofs.push_back(Mesh<2>::dofsOnGrid[k][6]);}
                    for(int k = seg[0] - 1; k < seg[0] * seg[1]; k = k + seg[0]){
                        Mesh<2>::boundaryDofs.push_back(Mesh<2>::dofsOnGrid[k][2]);
                        Mesh<2>::boundaryDofs.push_back(Mesh<2>::dofsOnGrid[k][3]);}
                    for(int k = seg[0] * (seg[1] - 1); k < seg[0] * seg[1]; k++){
                        Mesh<2>::boundaryDofs.push_back(Mesh<2>::dofsOnGrid[k][5]);
                        Mesh<2>::boundaryDofs.push_back(Mesh<2>::dofsOnGrid[k][4]);}

        
    return Mesh<2>::boundaryDofs;
};

//假定所有单元类型、阶数都相同

void Q2Mesh::distributeDofs()
{
    //以下示范DIM=2，_order_of_element = 1 | 2;
    //assert(DIM ==2);
    int k = 0;
    vector<Dofs<2>> _dofs;

    // for (Grid<2> &grid : Mesh<2>::gridList)
    // {
	    
	// _dofs.clear();
	// vector<Point<2>> p = grid.GetVertices();
	// _dofs.push_back({{p[0][0], p[0][1]}, k + k / seg[0]});
	// _dofs.push_back({{p[1][0], p[1][1]}, k + k / seg[0] + 1});
	// _dofs.push_back({{p[2][0], p[2][1]}, k + k / seg[0] + seg[0] + 2});
	// _dofs.push_back({{p[3][0], p[3][1]}, k + k / seg[0] + seg[0] + 1});
	// Mesh<2>::dofsOnGrid.push_back(_dofs);
	// k++;	    
    // }

    	for (Grid<2> &grid : Mesh<2>::gridList)
        {
            _dofs.clear();
            int a = seg[0];
            vector<Point<2>> p = grid.GetVertices();
            _dofs.push_back({{p[0][0], p[0][1]}, k / a * (4 * a + 2) + k % a * 2});
            _dofs.push_back({{(p[0][0] + p[1][0])/2.0, (p[0][1] + p[1][1])/2.0}, k / a * (4 * a + 2) + k % a * 2 + 1});
            _dofs.push_back({{p[1][0], p[1][1]}, k / a * (4 * a + 2) + k % a * 2 + 2});
            _dofs.push_back({{(p[1][0]+p[2][0])/2.0,(p[1][1] + p[2][1])/2.0}, k / a * (4 * a + 2) + 2 * a + 1 + k % a * 2 + 2});
            _dofs.push_back({{p[2][0],p[2][1]}, (k / a + 1) * (4 * a + 2) + k % a * 2 + 2});
            _dofs.push_back({{(p[2][0]+p[3][0])/2.0,(p[2][1]+p[3][1])/2.0}, (k / a + 1) * (4 * a + 2) + k % a * 2 + 1});
            _dofs.push_back({{p[3][0],p[3][1]}, (k / a + 1) * (4 * a + 2) + k % a * 2});
            _dofs.push_back({{(p[3][0]+p[0][0])/2.0,(p[3][1]+p[0][1])/2.0}, k / a * (4 * a + 2) + 2 * a + 1 + k % a * 2});
            _dofs.push_back({{(p[0][0]+p[2][0])/2.0,(p[0][1]+p[2][1])/2.0}, k / a * (4 * a + 2) + 2 * a + 1 + k % a * 2 + 1});
    	    Mesh<2>::dofsOnGrid.push_back(_dofs);
            k++;
        }

    // else
    // {
    //     std::cout <<"抱歉，干不了。\n"; 
    // }
};

int Q2Mesh::getTotalNumDofs()
{
    //此为Q_1_E算法，Q_2_E算法
    //return (seg[0] + 1) * (seg[1] + 1);
    // else if(_order_of_element ==2)
    	return (2 * seg[0] + 1) * (2 * seg[1] + 1);
    // else
    // 	return 0;
};


template<unsigned int DIM>
class UnstructedMesh: public Mesh<DIM>
{
public:
    // virtual void readMesh();
    // virtual void writeMesh();
    // virtual void distributeDofs(int _order_of_element);
    virtual void InputDofData() = 0;
    virtual void InputEdgeData() = 0;
    virtual void InputEleData() = 0;
};

class EasyMesh :public  UnstructedMesh<2>
{
private:
    string _MeshName;
    long int _ndofs;
    long int _nele;
    long int _nbndedge;
    vector<Dofs<2>> _DofsList;
    vector<Dofs<2>> _BndDofs;
    vector<vector<int>> _ElementData;
    vector<vector<int>> _BndEdges; 

public:
    EasyMesh();
    EasyMesh(string meshname);
    void InputDofData();
    void InputEdgeData();
    void InputEleData();
    int getNumGrid(){return _nele;};
    int getTotalNumDofs(){return _ndofs;};
    int getBndedgeNum(){return _nbndedge;};
    Dofs<2> getDofs(int i){return _DofsList[i];};
    vector<int> getElement(int i){return _ElementData[i];};
    vector<Dofs<2>>& getBoundary(){return _BndDofs;};
    //bool IsBoundary(int index){return !(std::find(_BndDofs.begin(),_BndDofs.end(),index) == _BndDofs.end());};
    vector<int> getBndEdge(int i){return _BndEdges[i];}; 
};

EasyMesh::EasyMesh(string meshname)
{
    _MeshName = meshname;
    InputDofData();
    InputEdgeData();
    InputEleData();
}

void EasyMesh::InputDofData()
{
    ifstream fin;
    fin.open(_MeshName+".n");
    string line;
    stringstream ss;
    getline(fin, line);
    ss.clear();
    ss.str(line);
    ss >> _ndofs;
    _DofsList.resize(_ndofs);

    int a, bndflag;
    double x, y;
    for (int i = 0; i < _ndofs; i++)
    {
	getline(fin, line);
	ss.clear();
	ss.str(line);
	ss >> a >> x >> y >> bndflag;

	Point<2> temp_point({x, y});
	Dofs<2> temp_dof(temp_point, i);
	_DofsList[i] = temp_dof;
	if (bndflag)
	    _BndDofs.push_back(temp_dof);	    
    }
}

void EasyMesh::InputEdgeData()
{
    ifstream fin;
    fin.open(_MeshName+".s");
    string line;
    stringstream ss;
    getline(fin, line);
    ss.clear();
    ss.str(line);
    int _nedge;
    ss >> _nedge;
    int a, d0, d1, bndflag;
    for (int i = 0; i < _nedge; i++)
    {
	getline(fin, line);
	ss.clear();
	ss.str(line);
	ss >> a >> d0 >> d1 >> a >> a >> bndflag;
	if (bndflag)
	{
	    vector<int> temp_edge(2);
	    temp_edge[0] = d0;
	    temp_edge[1] = d1;
	    _BndEdges.push_back(temp_edge);
 	}
    }
    _nbndedge = _BndEdges.size();
}

void EasyMesh::InputEleData()
{
    ifstream fin;
    fin.open(_MeshName+".e");
    string line;
    stringstream ss;
    getline(fin, line);
    ss.clear();
    ss.str(line);
    ss >> _nele;
    _ElementData.resize(_nele);

    int a, d0, d1 ,d2;
    for (int i = 0; i < _nele; i++)
    {
	getline(fin, line);
	ss.clear();
	ss.str(line);
	ss >> a >> d0 >> d1 >> d2;
	_ElementData[i].resize(3);
	_ElementData[i][0] = d0;
	_ElementData[i][1] = d1;
	_ElementData[i][2] = d2;
    }
}


#else
#endif
